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ABSTRACT 



A computer simulation study of channeling in a diamond 
lattice. The simulation was done for a xenon ion striking 
the (100), (110) or (111) surface of a silicon target. 

Potential functions for the Si-Si lattice bond and the Xe-Si 
interaction are postulated. The electronic stopping cross 
section for the (llO) channel of silicon is estimated. 

This work is a continuation in the development of a 
computer model formulated at the USNPGS which takes into 
consideration the displacement of the atoms in the target 
lattice as well as inelastic energy losses by the primary 
ion. The lattice was not thermalized and only the repulsive 
portion of the lattice-lattice potential was utilized. 
Computed ranges are in good agreement with experimental date. 
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INTRODUCTION 



I . 

Experimental and theoretical studies over the past decade 
have established that the slowing down of ions in crystals is 
strongly dependent on the crystallographic orientation of the 
lattice with respect to the direction of the incident beam. 

Ll, 2, 3,4,5] This effect is known as channeling. A particle 
moving along one of the major axis of a crystal may be "steered** 
(i.e. its momentum redirected toward the center of the channel) 
by successive gentle collisions with the **strings** of atoms that 
comprise the **walls’* of the channel. Davies et al [ l] reported 

125 

that this channeling effect increases the mean range of Xe 
by as much as a factor of 10 over the values observed in a 
non-crystalline solid of similar atomic number. 

Earlier computer studies performed at the NPGS investigated 
channeling in the body centered [6] and face centered cubic 
lattice. C 7 ] This thesis is an extension of this previous work 
into the diamond lattice. 

When a substrate is bombarded by a beam of incident ener- 
getic particles it will lose some of its own ions by sputtering 
and retain some of the incident ions. The incident ions retained 
are said to be implanted, and the technique of using an energetic 
ion beam to introduce ions into a substrate is called ion implan- 
tation. [5] This doping technique in silicon should be very 
effective due to the extremely open nature of its lattice and 
its large channel sizes. 
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The difference in the depth of penetration of energetic 
particles along the principal channels in silicon may be attri- 
buted to the various channel sizes. Figures 1-3 show the loca- 
tions of the various channels in the (100), (HO) and (111) 
surfaces of silicon while Fig. 4 is a "ball and stick" model 
showing the relative shape and sizes of the three principal 
channels in silicon. According to Davies [ l] and Manchester [8] 
the preferential order of channeling would occur along the ( 110) , 
(ill) and ( 100) channels in decreasing order with the (111) 
channel only marginally better than the (lOO) channel. This fact 
is graphically indicated in Fig. 5- 

Much work over the past several years has been devoted to 
studying this channeling effect in silicon. Manchester [8] has 
studied the feasibility of forming metallurgical junctions in 
silicon by making use of the channeling phenomenon. Eisen L 93 
has studied the channeling of medium mass ions through silicon 
with special attention directed toward the electronic stopping 
cross section of well channeled ions. Davies and his co-workers 
have written a series of papers [ 2 , 3 , 4>103 discussing the pene- 
tration of KeV and MeV projectiles in silicon and the disordering 
of the lattice caused by this ion bombardment, Glotin [ll], 
Nelson and Mazey [ 12 '] , and Dearnaley et al [ 13] have investigated 
the implantation of phosphorous in silicon. Eriksson et al C i4] 
have studied the implantation and annealing behavior of Group III 
and V ions in silicon. Gibbons f 5] sums up the current theoreti- 
cal and experimental work pertinent to the problem of predicting 
impurity distribution profiles in implanted material. 
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STUDY OBJECTIVES 



II . 

This report is an extension of the NPGS channeling model from 
the bcc and fee orientations into the diamond lattice. This 
thesis attempts to establish, with xenon as the bullet and sili- 
con as the target, the following characteristics for the three 
major channeling orientations: 

a) The potential between the bullet and the lattice, 

b) A value for the electronic stopping power, and 

c) The ion implantation profile at 5 KeV . 
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Ill . 



SIMULATION MODEL 



The basic model has been developed by the NPGS group over 
the past six years. A brief description of this model and the 
alterations made for this study follow but for detailed pre- 
sentations the reader is referred to Refs. 7, 15 and l6 . 

The simulation model consists of a single primary and a 
silicon (diamond) lattice target containing between 60 and 150 
atoms. All runs were made with the primary (bullet) striking 
the target normally, although the program is sufficiently 
flexible that the bullet could be fired into the target at any 
angle . 

The silicon target lattice as used in this simulation is 
composed of a central core and a surrounding shell. The program 
is structured so that the bullet is never allowed to leave this 
central core. In addition, in the (110) and (111) orientations 
the shell atoms were not permitted to move. This allows the 
program to run slightly faster and also tends to hold the core 
atoms in place for a slightly longer time. The shell atoms were 
not allowed to displace in the (100) channeling direction. Even 
with the shell atoms displacing, given a bullet energy and an 
equivalent channeling location in the (111) and (100) orien- 
tations, the (100) channeling program ran faster. For a detailed 
discussion of the lattice generator to include the dividing of 
the target crystal into active core and passive shell see 
Appendix A. 
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In previous computer simulations the interatomic potential 
function and force function were of the exponential (Born-Mayer) 
type 

F = exp (A + Bx) 

where A and B are empirically determined constants and x is the 
atomic separation. The functions are not always of this type 
in this simulation, and force and potential tables were con- 
structed. For details of this scheme see Appendix B. 

Two major assumptions are made in this simulation model, 

(1) The bullet moves in a perfect lattice undisturbed by the 
thermal displacements of the atoms. 

(2) Only repulsive forces between the lattice atoms have 
been included. 

While both of these assumptions are obvi'^usly erroneous, 
previous studies have shown that this simulation model is still 
a good first order approximation for channeling. 
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IV. PROCEDURE 



This simulation study was conducted in an IBM 36 O /67 computer 
using FORTRAN IV language. 

The computer sets up the diamond lattice in the desired 
orientation (see Figs. 1-4) and starts a xenon ion into it. A 
brief study of the symmetry of the target lattice shows that the 
indicated impact areas cover all the possible points for the 
three orientations of silicon. Within these impact areas, 
several points, depending on the orientation of the lattice, were 
selected as representative of the channel. In the determination 
of how many points were to be chosen in each impact area, con- 
sideration was given to obvious lattice symmetry and to computer 
running time. 

Various combinations of Bullet-Target Potentials and CELS, 
a constant used to determine inelastic energy losses, were tried 
at different bullet energies and the maximum range of the bullet 
for each orientation was matched against experimental range data 
as reported by Davies et al [ l] . 

Channel profiles for 5 KeV xenon atoms were determined for 
the (100), (110), and (111) channels of silicon. Runs were made 

to determine dE/dx for each major axis in silicon and maximum 
ranges for xenon at higher energies were determined knowing 




dE 



0 
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V, POTENTIAL FUNCTIONS AND INELASTIC LOSS CONSTANT 



The purpose of this section is to discuss first, the 
determination of the potentials used to approximate the Si-Si 
lattice bond and the xenon-silicon bullet-lattice interaction, 
and second, the determination of an inelastic loss constant. 
Presently, there exists in the literature calculated potentials 
for the Si-Si system [ I7l but none for the Xe-Si system, 

Eisen t 9] has published experimental electronic stopping cross 
sections for medium mass ions in silicon but no data exists for 
xenon in silicon. 

A, Si-Si BOND 

The starting point for selecting a potential for the Si-Si 
bond are the Interaction Potential Functions (IPF) as calculated 
by Harrison [ 17] . These potentials are only the repulsive portion 
of the potential function for pairs of particles and neglect any 
lattice effects. Harrison [ l7l also indicates that they are most 
likely too ”hard^^. Figure 6 shows the radial electron density 
for Si. Since this radial electron density was determined using 
a spherical approximation, it is not correct for r ^ 1,5a 

H 

(a = Bohr Radius = O. 529 S) [ 17] • Figure 7 shows several dif- 
H 

ferent potential functions for the Si-Si bond with the Gibson II 

for comparison. The Si-Si and Si-Si"^ IPF's show approximately 

the correct behavior and zero close to the nearest neighbor dis- 

+ 1 +1 

tance for silicon. The Si -Si also shows approximately the 
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correct behavior but fails to zero at the nearest neighbor 

^ . .+1 .+1 

distance. Curve 3 of Fig- 7 is the Si -Si IFF forced to 
zero at the nearest neighbor distance. 

Any one of these three potentials (Si-Si, Si-Si or Si-Si 
(zeroed)) will serve as a good approximation for the Si-Si bond. 
It is obvious that the Gibson II will not serve as a good repre- 
sentation of this bond. 

The covalent bonding between Silicon atoms has the effect 
of concentrating about 0.5 of an electron of the outer shell 
along every chemical bond Cl8l. Therefore, it seems reasonable 
to conclude that each silicon atom along the Si-Si bond sees the 
other in a slightly ionized state but not completely ionized 
to the +1 state due to the shielding effects of this shared 
electron. This would seem to eliminate the Si-Si and Si"*"^ -Si"*"^ 
potentials for the Si-Si bond and lead to the selection of the 
Si-Si"*" potential as representing the bond. It was noted earlier, 
however, that the IPF^s are probably too hard. Therefore, the 

Si"*"^-Si'*'^ (zeroed) potential was chosen as representative of the 

. . + 

bond potential as it is "softer” than the Si-Si potential and 
is forced to zero at the nearest neighbor distance. 

Admittedly the reasons for choosing the Si -Si potential 
were rather arbitrary but the reader is asked to keep three 
points in mind. First, all three potentials are approximately 
equal and trial runs indicated no significant effect on the range 
when these potentials were interchanged- Second, the neutral 
silicon atoms are too large to fit into the Si lattice. Third, 
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the two overwhelming factors that control channeling ranges are 
the bullet-lattice potential and the inelastic loss approximation. 
Any reasonably potential selected for the Si-Si bond should have 
no discernable effect on these ranges. 

B. INELASTIC LOSS CONSTANT 

One of the principal mechanisms of energy loss for channeled 
ions in a solid is the interaction of the ion with electrons 
(both bound and free) of the solid. This computer study accounts 
for this loss through the use of a "frictional force multiplier" 
which acts on the ions velocity in a manner to produce a "drag" 
force on the incident ion. This "friction force multiplier" is 
directly proportional to the electronic stopping cross section. 

Eisen [9] reports electronic stopping cross sections (S^) 
for medium mass ions (5 ^ Z ^ 19) in the <100), <110), and <lll) 
channels of silicon. In unpublished data [ 19] he extends for 
the (no) channel by allowing the atomic number of the incident 
ion to increase up to Z = 36 (Krypton). Eisen’s data reveals the 
following: first, was found to exhibit a strong oscillatory 
dependence on Z, the atomic number of the incident ion, second, 
is lowest in the (H0) channel and third, is approximately 
equal in the <lll) and (100) channels. 

No data presently exists in the literature for the electronic 
stopping cross section of xenon in silicon. However, based on 
previous computer studies of channeling in copper and tungsten 
undertaken at NPGS the frictional force multiplier was estimated 
to be 1.0 X 10 n-sec/m in the { lOO) channel. Using Eisen 's 
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data [ 19 ] and noting the ratio between in the (100) and (110) 

channels, the frictional force multiplier in the (llO) channel 
is approximately -8 x 10“^^ n-sec/m. These values were used in 
this computer study, 

Ericksson et al U^O] give estimates of electronic stopping 

cross section for various ions along the (100) and (llO) channels 

of tungsten. Their data also indicates that S shows a strong 

e 

oscillatory dependence on Z, the atomic number of the incident 
ion and. -this oscillatory dependence is very similar to the 
dependence Eisen [9?19] found for ions channeled in silicon. If 
the assumption is made that the ratio between stopping cross 
sections for xenon and krypton (which is chemically similar to 
xenon) channeled in tungsten is approximately the same as the 
ratio between these ions channeled in silicon, one obtains a 
value of 

2 

8.5 10 ^y . (where v = 1.0 x 10^ cm/sec) 

atom ' ’ 

as an estimate of S in the (110) channel in silicon. This 

e 

“i4 

converts to a value of 7 10 n-sec/m as a frictional force 

multiplier for the (llO) channel and is in good agreement with 

the value used in this computer study. 

^ ^e (v)*TI*l-6 X 10 j/eV) , ^ ^ 

CELS = - ^ ^ ^ ^ ^ where T] = number density. 

v( cm/sec) x 10 (m /cm ) 

The author does not make any claim that the values of the 
’’frictional force” used in this study are ’’correct”. It is felt, 
however, that based on a study of Eisen *s data [9,19] smd the 
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extrapolation described above, the values chosen for S will 

e 

certainly be within 20% of the experimentally determined values 
and more than likely will be much closer than that. 



C. Xe-Si POTENTIAL 

As in the case of the Si-Si bond, no potential function 
exists for the xenon-silicon interaction. Therefore, the 
Interaction Potential Functions will again serve as a starting 
point for selecting the Xe-Si potential. Figure 8 shows the 

radial electron density for xenon in the nonionized and +6 

. . . . . +6 .+4 

ionization state, While Fig. 9 shows the IPF for Xe -Si and 

+6 . +1 
Xe -Si . 



Far from the Si-Si bond, the electron density for silicon 
should approach the Si"^^ distribution (Fig. 6). Special note 
should be taken of the (^111) channel of silicon where the atoms 
bounding the channel have one of their bonds parallel to the 
centerline of the channel. Determing the ionization state of 
xenon is also difficult. Harrison [ 17] has found for Xe-W 
system a good potential fit can be found if xenon is assumed to 

be in the +6 or +8 ionization state, even though this high ionic 

+ +6 . ^ , 
state IS somewhat surprising. The Xe -W potential was found 

to be unrealistically strong. As a result, of the above it was 

+6+4 . . 

concluded that the Xe -Si IPF was a good starting point for 

the bull et -la ttice potential. 

The Xe^^-Si^^ IPF (Fig. 9) diverges to unrealistically 
strong values at large separations. This divergence is charac- 
teristic of the theoretical model used to calculate the IPF 
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and is not characteristic of the physical system [ 17] • A Born- 
Mayer approximation for the Xe-Si interaction would fail at small 
separations because it would approach a finite value at zero 
separation, while the actual Xe-Si function should approach a 
pure coulomb potential at small separations. Our final choice 
of the potential to represent the Xe-Si system incorporated both 
of the above approximations. For separations less than 1.2a 

H 

the IFF was used, while a Born-Mayer approximation was constructed 
for separations greater than 1.2a . 

Curve BMl (Fig. 9) when used in conjunction with the pre- 
viously discussed values of the ^’frictional force" parameter 
gives excellent results for xenon ions injected in the ( 110) and 
(100) channels, however, it is a poor approximation for xenon 
in the (ill) channel. After several trial runs it became obvious 
that all three channeling directions could not fit on one Born- 
Mayer approximation. The maximum range in the (ill) channel 
could be forced to approach the experimental range if the 
"frictional force" constant is decreased by two orders of magni- 
tude. This is unrealistic since it conflicts dr ast ically w ith 
Eisen’s [9] experimental data. 

A solution to the problem was to make a separate Born-Mayer 
approximation for the (ill) channel; this is curve BM2 in Fig. 9* 
Intuitively, one feels that this answer is not the correct one. 

All previous studies have found that one potential satisfies all 
channeling directions. However, two points should be kept in 
mind. First, all previous computer studies dealt with channeling 
in a metallic target where the principal method of bonding between 
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lattice atoms is metallic. This is the first computer study where 
the bonding between the lattice atoms is co-valent and any con- 
clusions drawn for a metallic target need not necessarily be true 
in silicon. Secondly, the orientation of the bond between silicon 
atoms with relation to the channeling directions varies among all 
three principal channels of silicon. As previously mentioned, 
the extreme is reached in the (ill) channel where the Si-Si bond 
is parallel to the center line of the channel. This more than 
likely has an effect on the silicon charge distribution as seen 
by the xenon ion and could possibly cause the potential for the 
<111) channel to differ from the other channels. 

Another possibility is that a Born-Mayer approximation is not 
valid in silicon. A picture of the Xe -Si potential might be 
obtained by combining the portion of the three curves: IFF, BMl 
and BM2 that we know give us accurate results. Such a composite 
curve is shown in Fig. 10. It is obvious that this composite curve 
is unrealistic. 
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VI . RESULTS 



A. <100> CHANNEL 

eV 

Figure 11 is a plot of dE/dx( in /LU) versus ion energy for 
a perfectly channeled xenon ion in the (100) direction. (Note 
ILU = 2 . 71 ^ = 0 ,0646 ^g/cm^ ) . Between 3 O-IOO KeV, dE/dx is 
approximately proportional to the energy indicating that in this 
region the elastic and inelastic collision processes are both 
exerting a significant influence. This agrees with experimental 
data reported by Davies [ l] . 

Figure 12 is a plot of dE/dx versus E . Above oOO KeV 
dE/dx IS proportional to E indicating that as expected the in- 
elastic process is the primary loss mechanism. Both Figs. 11 
and 12 show that below 20 KeV the elastic collision process is 
primarily responsible for slowing down the ion. 

Figures 13 and l4 are range profiles and range profile 
sections for the (100) channel. Figure 13 shows only one-half 
of the channel, the range profile for the other half is very 
similar but not identical to this. 

Table 1 (below) is a comparison between experimental ranges 

obtained by Davies [l] and computer calculated ranges for the 

( 100 ) channel. The column labeled dE/dx is the range value 

obtained by using the dE/dx curve in Fig. 11 and the relationship 

E 

R =j ( dE/dx) dE . 

0 
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The eomputer ealeulated range at ^0 KeV, 


although higher than 


the 0 cper imental 


value reported 


by Davies 


, falls within his ^3% 


error estimate. 










Range ( lOO) 




Energy 

5 


Exper imental 


dE/dx 

60 


Computer 

69.7 


20 


491 . 7 * 


485 


490,2 


40 


1162.7 


1230 


1193.8 



Table !• Energy in KeV, range in LU 

*Exper imental range for 20 KeV extrapolated 
from Davies [4] experimental data. 

In a bcc or fee erystal an ion injeeted into the eenter of 
the ehannel remains there • This is not the ease for an ion 
injeeted into the eenter of the (100> ehannel of silieon. The 
ion "wobbles** slightly around the ehannel eenter oseillating 
between (110) planes. (Fig. 15) Note the similarity of the 
orbit in the seeond and fourth quadrant. 

Figure 6 is a Range Profile Seetion for a 5 KeV silieon ion 
injeeted into the (100) ehannel of silieon. Although no experi- 
mental data exists to eompare maximum ranges it is interesting to 
note that the range profile seetion is very similar to the range 
profile of xenon in silieon. 

B. (110) CHANNEL 

Figure 7 is a plot of dE/dx versus energy at various values 
of CELS (analogous to eleetronie stopping power) for a xenon ion 
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injected into the center of the (llO> channel. It is apparent 
that the energy loss is very sensitive to CELS especially at 
higher energies indicating that electronic stopping is the 
principle mechanism of energy loss for well channeled ions in 
the (110) direction. Figure l8 is a plot of dE/dx versus 
(ion energy) for well channeled ions with CELS = - 8 x 10 
n-sec/m, the value selected for use in this study. dE/dx is 
proportional to E at incident ion energies greater than 2 KeV 
again indicating that electronic stopping is the principle loss 
mechanism . 

Figure 19 and Fig. 20 are a range profile and a range profile 
section for the (110) channel. Both show only half the channel 
although the other half of the channel should be identical. 

While it is apparent that the (110) channel is a preferred 
channeling direction the maximum range falls off rapidly if the 
ion is injected slightly off center. (For example, if the ion 
is injected approximately 0.2LU from the center of the channel 
the maximum range decreased by one^-half ) . Of particular interest 
is point 1 in the Range Profile Section (Fig. 20). An ion in- 
jected at this point in the channel has a range of lOOLU which 
is much greater than expected. The xenon ion injected here is 
strongly influenced by the silicon atom located at position C 
which forces it into the center of the channel and adds signi- 
ficantly to its range. 

Table II is again a comparison between ranges obtained 
experimentally by Davies [ l] and computer calculated ranges. 
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Range 110 



Ener gy 


Exper imental 


dE/dx 


Computer 


5 


— 


1090 


1094 


20 


2371 


2400 


-- 


40 


3812 


3800 






Table II. Ranges in LU, Energies in KeV. 

The ranges calculated using dE/dx are in good agreement with 
experimental ranges obtained by Davies [ l] • No attempt was made 
to channel a 20 KeV xenon ion as it would be too expensive in 
computer running time (approximately 90 minutes). 

Figure 21 shows, the ranges obtained for selected points 
in the ( 110) channel when the energy of the incident xenon ion 
was increased to 20 KeV. Two points should be noted. First, 
the channel ^expands” at higher energy as expected. Second, a 
comparison of the range of the ion incident on point 1 with the 
range obtained when the incident ion energy was 5 KeV (Fig. 20 ) 
shows that the large range for the 5 KeV ion is a low energy 
phenomenon . 

C. < 111 > CHANNEL 

Figure 22 is an enlargement of the (lH) channel. Initially, 
point 1 was thought to be the center of the channel but subse- 
quent probing of the channel proved this assumption to be incor- 
rect. The points numbered 2 and 3 3 iXe actually the channel 
"centers". If atoms A, B, D and B, C, D are thought to form 
triangles, as indicated, then points 2 and 3 t^e considered 
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the center of their respective triangle equidistant from the atoms 
that form the triangle’s vertices. Additionally, although the 
basic shape of the channel remains unchanged, the atoms that 
bound the channel do not necessarily lie in the planes indicated 
in Fig. 22. Two other possible configurations ares (a) atoms A 
and C in the front plane, B in the third plane and D in the fifth 
plane or (b) atoms A and C in the fifth plane back, D in the third 
plane back and B in the front plane. These three channel con- 
figurations are shown in Fig. 3. They affect the ion only as it 
initially enters the lattice, by determining in which direction 
it will be initially steered and have almost no effect 0.5%) 
on the calculated maximum range. 

Figure 23 and 24 are plots of dE/dx vs ion energy, E and 

... 

E . Up to 100 KeV both elastic and inelastic collisions are 
both important loss mechanisms. Below 20 KeV, elastic losses 
appear to predominate. 

Figure 25 shows dx/dE curves for xenon in the (111} channel. 
The area under these curves should be an accurate representation 
of the total range travelled by an incident xenon ion. The solid 
curve (I) was obtained by injecting xenon ions into the center of 
the (111) channel at various energies, allowing them to penetrate 
the lattice for a short distance 5LU) and determining the 
energy loss of the xenon ion per lattice unit. The dotted curve 
(II) is actually a 20 KeV run for a xenon ion initially injected 
in the center of the (111) channel and followed until it finally 
stops in the crystal. The difference between these two curves 
can be accounted for by noting that an ion injected into the 
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center of the (111) channel deviates considerably from the 
channel center (up to .09LU) as it passes thru the crystal and 
that this deviation increases as the ion energy decreases. 
Therefore, the implicit assumption made in determining curve I, 
that the ion remains in the center of the (ill) channel as it 
passes through the crystal is incorrect and will introduce large 
errors (up to 20%) in ranges at the low end of the energy spec- 
trum where the elastic loss mechanism predominates. Table III 
is a summary of maximum ranges obtained at various energies for 
the (ill) channel. 

Range (ill) 



En er gy 


Exper imental 


Curve I 


Curve II 


Computer 


5 


-- 


37 


60 


57 


20 


511 


417 


500 


511 


40 


1209 


1050 


1170 


^ - 



Table III. Energy in KeV, range in LU. 

Figure 26 is a representation of the (ill) channel and shows 
ranges obtained at different initial impact points in the channel* 
Two things can be noted about the range distribution. First, 
there is a definite peak at the center of the channel and 
secondly, the total range is extremely dependent on exact im- 
pact point in the channel. This last fact makes it impossible 
to draw any meaningful range profiles or ranges profile section. 
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VII « CONCLUSIONS 



1. Maximum channeled ranges were fit to experimental data 
at 20 KeV, 

2. For the ( 100) and (ill) channels both elastic and 
inelastic losses are important between 20 and 80 KeV. 

3. Inelastic losses become predominant at energies greater 
than 400 KeV for the ( lOO) channel. 

4- Inelastic loss seems to be the predominant loss mechanism 
for well channeled ions in the (110) channel even at low energies 
(--2 KeV). 

5» The range in the (ill) channel is extremely sensitive to 
impact parameter and appears to possess no symmetry. 

2 

6, S was determined to be 8.5 x 10 ^ eV cm (llO) 

e atom ' 

channel of silicon. In the ( lOO) and (ill) channels 

2 

S = 10.6 X 10 eV cm , _ ^ ^^8 

e atom ^ 

7c All three channeling directions for xenon in silicon can 
not be determined using the same Born-Mayer Potential Function. 

A minimum of two Born-Mayer Potentials are required. 

8. At low energies the (100) channel is the preferential 
direction for channeling when compared to the (ill) channel. 

This is the reverse of the situation for incident ions energies 
greater than 20 KeV. 
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9. Although the potentials and the ’’frictional force” con- 
stant were discussed separately it can not be overemphasized 
that all three of these variables must be treated as an inte- 
grated whole. A change in either the bullet-target potential 
or the frictional force constant has a profound effect on the 
maximum channeled range especially in the ^ llO) orientations. 

The final selection of the Born-Mayer approximations to the 
Xe-Si interactions, the frictional force parameter and Si-Si 
bond potential seems to represent the best of all the parameters. 
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APPENDIX A (LATTICE GENERATOR) 



Previous computer studies undertaken by the NPGS Group over 
the past several years have dealt exclusively with the body 
centered and face centered cubic lattices. In building the 
diamond lattice, which is basically two interpenetrating face 
centered cubic lattices, the primary aim was to construct the 
lattice in such a way that only mimimum modifications on the 
basic computer model for channeling would be necessary. The 
diamond lattice generator was constructed so that no major change 
in the basic lattice regenerator was necessary. Therefore, to 
understand the lattice generator, it is necessary to understand 
the fundamental mechanism behind the regenerator subroutine. 

Subroutine REGEN works as follows: 

(1) It receives the parameter MSS from the main program 
which indicates the direction of regeneration (_^X, j^Y, _+Z) . 

(2) It checks every atom in the crystal against a constant 
which depends on the direction of regeneration and either 

a, discards the atom since the bullet has passed and is no 
longer under its influence, or 

b. shifts and renumbers the atom while preserving its 
velocity components. 

(3) It builds and undisturbed section crystal in the 
direction of travel of the bullet. 

( 4 ) It shifts the bullet to its proper location in the 
rebuilt crystal preserving its velocity components. 
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Schematically, this is done as follows: (The crystal regen- 



erates when the bullet passes plane E) 
A) Crystal Prior to Regeneration 
Direction of Travel 




Disordered Crystal, 
Passed By Bullet 



B) Crystal During Regeneration 




The foilwing planes 
are identical: 

A - E - I 
B - F 
C - G 
D - H 



Plane E goes to A 

F »» »» B 

G C 

H ” '' D 

I ft ft E 

Old planes A,B,C,D, 
discar ded 



Discar ded 



C) Crystal After Regeneration 



New 

Plane 



Was 

Plane 






-bullet 



B C 


: I 


-) 1 


1 


' c 


; H 


F C 


5 1 











Planes F,G,H,I are newly 
constructed as part of a 
perfect crystal lattice. 



This scheme discards the portion of the crystal the bullet 
,has passed through and builds an undisturbed portion of the 
crystal in front of the bullet. 



31 



To build a crystal which will fit into the regeneration 
scheme the crystal generator must: 

(1) Construct one plane at a time in the Y direction. It 
places the first atom in the Y = 0 plane at the lowest value of 
X and Z (0,0) and continues placing atoms at increasing values 
of X for fixed Z? then increases the value of Z and repeats the 
process until the plane is completed. The next Y plane is con- 
structed in a similar manner and additional Y planes are 
constructed until the microcrystallite is completed. It assigns 
a number to each atom. A completed Y = 0 plane for the (100) 
orientation has atoms positioned and numbered as follows 




Planes must be generated in this manner. 

(2) Build the following symmetry into the crystal. The left 
and right face are identical to each other and identical to one 
plane in between them. There is similar symmetry between the 
front and back faces as well as the bottom and top faces. 

This symmetry in the lattice allows the regenerator to 
operate as follows: 

(1) It checks the location of each atom in the crystal to 
see if it is to be regenerated (shifted). 
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(2) If it is to be regenerated, a fixed number (IDX, IDY, 

or IDZ depending on the direction of regeneration) is subtracted 
from the atom number. This partially shifts the atom into its 
identical plane. 

( 3 ) The regenerator then subtracts (or adds) DXT , DYT, or 
DZT from one of theatom’s co-ordinates depending on the direction 
of regeneration. This completes the movement of the atom into 
its identical plane. 

( 4 ) As previously mentioned the regenerator preserves the 
velocity components and potential energies of all atoms 

r egeher ated. 

( 5 ) The bullet is regenerated in a similar (but not identical) 
manner , 

In addition to being tailored to fit SUBROUTINE REGEN the 
lattice generator, through the use of a variable LCUT, breaks 
the lattice into an active core and a surrounding shell. This 
is an absolute necessity due to the large size of the crystal 
which was dictated by the symmetry requirements of SUBROUTINE 
REGEN. (Orientation (100) contains 96 atoms, orientation (110) 
contains 60 atoms and orientation ( 111 ) contain^s l 46 atoms). 

LCUT does this in the following way: 

1) Any atom in the active portion of the crystal is given 
an LCUT = 0, This allows it to interact with any other atom 
in the crystal . 

2) Any atom in shell is given an LCUT = 1 which allows it to 
interact only with core atoms and not with any shell atoms. 
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The lattice generator subroutine also contains the constants 

XLL, XLFj ZLLj ZLF which are used to determine when the lattice 

is to be regenerated. They are chosen to insure that the bullet 

always remains in the core and the lattice regenerates when the 

bullet is about to leave the core and enter the shell, 

125 

As an example, a 20 KeV Xe atom in the center of the (100> 
channel of silicon has a range of 491*7 bU and the program runs 
for 35 minutes and 39 secs. With the LCUT package inserted into 
the program the calculated range is 490.2 LU with a program 
running time of I 3 minutes and 12 seconds. The LCUT package is 
obviously effective for it conserves computer time without 
sacrificing any accuracy in range calculations. 

In summary, the Lattice Generator 

(1) Contains all the constants necessary to regenerate the 
lattice , 

(2) It constructed, with appropriate tests, to insure that 
it fits into the scheme necessary for Subroutine REGEN, and 

( 3 ) Contains the LCUT package which reduces the number of 
active atoms in the crystal, thus conserving running time. 

An alternative to building the lattice generator to fit 
subroutine REGEN would be to construct a smaller active crystal, 
eliminate the LCUT package and rewrite the REGEN subroutine to 
fit the lattice generator. This method was not chosen for the 
following reasons: 

(1) It was easier to modify the existing lattice generator 
to fit the old Subroutine REGEN than to start from scratch and 
build a new Regeneration Subroutine and Lattice Generator Sub- 
routine tailored for each other, and 
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(2) Simple modifications in the generator alone would not 
drastically alter the existing Channeling Program which had 
proved to be an accurate representation of the channeling pro- 
cess in fee and bcc crystals. 
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APPENDIX B (CALCULATION OF FORCES AND POTENTIALS) 



A major change in procedure was introduced in the calculation 
forces and potentials. In previous models the interatomic po- 
tential function and force functions were of the exponential type, 

F = exp (A + Bx) 

where A and B were empirically determined constants and x is 
the atomic separation. The atomic separation was determined, 
substituted into the appropriate function and the force and 
potential determined. 

This computer simulation determined forces and potentials as 

follows. First, a set of potentials (in Rydbergs) derived from 

Interaction Potential Functions are read into the program for 

the lattice-lattice and bullet -la tt ice interactions. These 

potentials start at 0,01a (a = Bohr radius = 0,529^) and are 

H H 

spaced every 0,01a up to 5«*0a interatomic separation. Secondly, 

H H 

a control card associated with each deck of potentials is read. 
This control card allows the extrapolation of a Born-Mayer type 
potential (if desired) from any point on the Interaction Potential 
Function Curve. Thirdly, the potential at the nearest neighbor 
distance is determined and then the potential is zeroed at that 
distance. All potentials at atomic separations smaller than the 
nearest neighbor distance are adjusted accordingly. All poten- 
tials at separations larger than the nearest neighbor distances 
were set equal to zero. This does not leave any potential free 
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regions in the lattice. Fourthly, potentials are converted to 
electron volts. Finally, the relation between the force and 
potential energy can be expressed simply as 

F(X) = - VU(X) . 

Assuming that the potential is linear between any two successive 
points in the potential table, the force can be expressed simply 
as 

Force = (U(X) - U(X + A)) -(const) 
where the constant converts the force to proper units. 

The determination of a force or potential from the tables is 
Just a simple matter of interpolation. For example, if the atomic 
separation between two atoms is 1.2752^^ would Just make a 

linear interpolation between the values at 1 , 27 a and 1 . 28 a . 

H H 

These tables are used in both SUBROUTINE STEP which calculates 
forces and SUBROUTINE lONPE which calculates potential between 
the ion and the atoms . 

Two points should be noted. First, the nearest neighbor 
separation is approximately 4*436a . Therefore, for a separation 
of 4-434a the interpolation is between the value in the table 
at 4*43a„ and the value at 4*436a . Secondly, the average force 
between two points is assumed to occur halfway between them. The 
interpolator for the forces accounts for this fact. 

One major problem remains unsolved in the construction of the 
force table (not in the interpolation of forces in SUBROUTINE STEP 
after the table has been constructed). While the potentials 
steadily increase as the interatomic distances decrease, the 
forces in the force table do not monot onically increase with 
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decreasing atomic separation. This is true only when potentials 
derived from the Interaction Potential Functions are used. No 
difficulty arises in the region where a Born-Mayer type of ex- 
ponential fit for the potentials is made. Here the forces behave 
as expected. 

The reason for this discrepancy is found in the calculation 
of the interaction potentials that are used by the computer pro- 
gram as data. The interaction potential is calculated exactly 

from 0.01 to 0.20a and every 0.20a to 5*0a (i.e. 0.2, 0.4, 

H H H 

0.6a would be calculated exactly, 0.50a would not). Potentials 

H H 

in between these values are calculated by linear interpolation. 
Hence, the forces in the table would remain essentially constant 
in the region where the potentials were derived by linear inter- 
polation and would only change abruptly when a point is reached 
where the potential was calculated exactly. Thus while the 
overall force would increase as interatomic separation decreased 
there would be regions where the force stayed essentially con- 
stant or even decreased with decreasing atomic separation. 

Two main efforts were made to eliminate this discrepancy both 
of which began by considering only the potentials that were cal- 
culated exactly and discarding those that were determined by 
linear interpolation. First, subroutine ALI and ATSE from the 
scientific subroutine package of IBM were used to fit known po- 
tentials with a polynomial and using this polynomial, points in 
between the known potentials were interpolated. In the second 
method an exponential fit of the known potentials was made, and 
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an interpolator constructed as follows. Assume the potential P 
at two points x and x are known and are to be joined by an 

-L ^ 

exponential. We desire to know the potential at x^ . 




2 1 

AX = — - — where N is an integer 

(1) = Aexp(-Bx^) 

(2) = Aexp( -B(x^+x^) ) or P^ = P^ exp(-Bx^) 

Similar ly 

(3) P 2 = P^exp(-BNAX) or ln(P 2 /Pj^) = -B*NAX) 

(4) BAX = - ln(P2/Pj^) )/N 

Assume X = M*X where M is a constant 0 then 

a 

= P^ exp(-BAX'M) 

Both the polynomial and exponential fit of known potential 

yielded approximately the same results. In general, for atomic 

separations of less than 1.6a the forces increase with decreasing 

H 
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atomic separation as desiredo For separations greater than 1.6a 

H 

the forces act as expected, except that in the immediate vicinity 



('^.02a ) of a calculated potential the forces decrease slightly 
and then increase again. Since there was little to choose between 
the two methods, the exponential interpolator was selected. The 
interpolator constructs the tables slightly faster (~ 2 secs) and 
requires less storage space. 

The effects of this problem on the force tables are as 
follows : 

+6 .+4 

a) Xe -Si -no effect. The forces behave as expected since 

the Born-Mayer approximation is constructed beginning at 1.2a . 

H 

- “^4 • "^4 .... 

b) Si -Si (used as bull et -latt ice interaction potential 

in preliminary studies of (100) channeling). No effect, as a 
Born-Mayer approximation is constructed beginning at 1.2a . 

c) Si'^^-Si'^^ (used as lattice-lattice interactions for all 

directions of channeling). This force table contains the error 

discussed above for interatomic separations between 1.4 and 

3»0a . However, this error has negligible effect on calculated 
H 

channeled ranges since (a) there is no error in the force table 

at atomic separations of 4»4a , the approximate nearest neighbor 

distance silicon atoms and (b) the bullet passes the lattice atoms 

so quickly that when two lattice atoms have an interatomic spacing 

of less than the bullet is no longer under their influence. 

H 

While this error in the force table has virtually no effect 
on channeling range estimates, before the force table concept is 
used for other computer studies this problem must be solved. 
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APPENDIX C (PROGRAM PARAMETERS) 



The following parameters have been added to the basic computer 
channeling model (Chan 69MODO). 

BPNT - Bullet-La tt ice Potential Energy 
TPNT - Lattice-Lattice Potential Energy 
BFORCE - Bullet -Lattice Force 
TFORCE - Lattice-Lattice Force 

KA - an integer constant which identifies the point in the 

potential table where a Born-Mayer potential of form, 

-Bx 

Ae , begins 

BPCNST - the value of B in the Born-Mayer Potential approximation 
for the Bullet -Lattice Potentials 
TPCNST - the value of B in the Born-Mayer Potential approximation 
for the Lattice-Lattice Potential 
ICHECK, ICOUNT, NCOUNT - constants used in the exponential inter- 
polator between calculated interactions potentials in 
constructing the Potential Tables 

BDELX - a variable (BAX) (see appendix II) used in the exponential 
interpolator between interaction potentials in con- 
structing the Potential Table 

ITITLE - 80 alpha numeric characters used on a heading card for 
potential data deck 

LCUT - integer constant with value of either 0 or 1 • Used to 
determine which atoms are in the core (0) or shell (1) 
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of the crystal 



CVEDA - conversion factor used to convert electron volts per 
Bohr to newtons per meter 

CLU - conversion factor used to convert distances in LU to 
distances in units of 100a 

H 

AFC - distance in 100a 

H 

IFC - integer value of AFC 

CLU, AFC, and IFC are used in the following way. Assume 
interatomic spacing between atoms = 0 . 5LU 
AFC = (CLU) X (DIST) = 256.14 
IFC = 256 

An AFC of 256,14 corresponds to 2.56l4a 

H 

An IFC of 256 corresponds to 2.56a 

H 

AFCl “ nearest neighbor distance in units of 100a 

H 

AFC2 - the reciprocal of the difference between the nearest 

neighbor distance and 443-0. Used in Subroutines STEP 
and lONPEo 

DFFl - the difference in 100a between the atomic separations of 
the bullet and lattice atom and the nearest neighbor 
distance of silicon 

GETIME, SETIME - machine subroutines used to determine actual 

program running time. They terminate program execution and 
have cards punched before exceeding the running time allo- 
cated to the program by the computer. This allows a re- 
start capability. 
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Computer Program - SUBROUTINE EX^lOO 



SUBROUTINF OLIOO ♦♦♦♦♦♦♦♦♦♦♦♦♦>)(»♦♦>)(♦♦♦ 

♦ THIS IS 4 LATTICE PENERATHR FOR TUP 01 AMONU Lft'^'^ICF IN ♦ 

* THE (IOC) fRIENTATIPN. * 

♦ '^HE CRYSTAL OpvFLHPS IN TP'S ORDER XfFOLLOWED BY 7, ♦ 

* PnLLPWFD BY Y * 

COMMON /CP^'1/RX( 200) , R Y ( 2 00 ) , R ? ( 200 ) ,LCUT (200 ) ,1 L,LD 
common /C r M2 /R OF ,RCF2 fROFM, AC , PAC, PPTC, FTC, DFP’T , F-’TO, 

2 FM,PPIV,TR0T 

CnMMPM/crv4/I X,I Y, I 7 , 1 XP.IYP, I7D,SCX ,SCY ,SC7, Fo 

COMMCN/CCM7/RI ,LSS,SPX ,SO7,C0X,C0Y ,C07 

common / ccMq/ i['x,iOY,in7 (200) ,nxT,nYT,rzT,TDcx,'rroY» 

2 TP07 

COMMON /CCM22 /XLL ,X(F,7LL,7LF 

« ^HF FOLLOWIKO CONSTANTS ARE USED FITHFR BY THIS SUP- * 

♦ PnuTTNF TO OFNERATF THE LATTICE OR IN OTHER SURRCHT INFS * 

* AS RFoiiIPEO. ♦ 

XLL=C.CQ 
XLF=3.Cl 
7LL=C. SP 
71 f = 3.01 



IB = l 
IDX=1 
IDY=41 
T X=F 
I Y=B 
17 = '^ 

PCX=1 .0 
SC Y= 1 . C 
SC7 = 1.0 
nXT=2. C*PC X 
DYT=2. C*SC Y 
DZT=2.0'('SC7 
DY=0. 5 
nx=o .5 
07 = 0.5 

* the following serifs of NESTEO OC loops r,FNERATE -^HF ♦ 

♦ LAT'riCE ANO INSURE IT CONTAINS THE REQUIRED SYmwc-^ry » 

JT=0 
Y=- sc Y 

DO 6C J=l,IY 
Y=Y4-SC Y 
KT=0 
7=- SCZ 

DO 50 K=1,I7 

z=z+scz 

IT=0 
X=- SC X 

DO 4C 1=1, IX 

X=X-*-SC X 

ITT = I T^-JT^-KT 

IF( ITT-( ITT72*2) ) 40,30,40 
30 RX(M) = X 
oy(M)=Y 
R7( M)=7 
ID7 ( M )=5 
L = M 

m = m+ 1 

40 I'^=IT+l 
50 KT=KT+1 

IF( J.E0.5) GO TO 80 
on 200 IA=LA,L 
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TP( A.C^SCX-RXnAJ^SCX) 200,200, ?0? 

202 TP( 7-RM lA) ♦SCZ) 200,200,203 

203 R X{ I +TP )-R X( lA ) +nX 
PY( U-TP )=RY( lA )+DY 
R7 ( L+IP )=P7( IA)+07 
in7( L+TP 

IB=IP+1 
200 C'^NTINIIF 
M = (_ + C 

LA 

TP=1 
JT=JT+1 
PC LL=M-1 

* the POLl.nwiNG or I OOP DIVI^fs the CRYSTAL INTO AN ACTIVE* 

* CORE AND AM INACTIVE SURROUNr'ING SH«=LL. * 

on 2 T=2,LL 
LCIIT( I »=0 

IP(P X( I ) .CT. :j. l.OR. Rxn ) . L^.O.^ ) LCUTd ) = 1 
TF(P Y( T ) .CT.2. 5) LCUT{T)=1 

IP(P7(n.l,F.C.9.'^R.R7{I).r,':.-^.i) lCUT( t ) = i 
2 CT'NTTNDF 
LCUT( 1 ) =0 

R 1=( Rpr+DY( 2 ) ) /APS (COY) 

PE'^IIPN 

FN.r 
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Computer Program - SUBROUTINE CLllO 



SUBPOUTINE r>LllO ♦♦♦«♦#♦*♦♦♦♦♦♦♦♦♦♦♦♦* 

♦ THIS !S t LATTICE f7ENPR/M0R PQR THE DI&MQND LATTICE IN ♦ 

♦ the ( 110 ) CRI'^NTATI ON. ♦ 

♦ the CPYST?l ncvELHPS IN THE CROER X.PCLLOWFr RY 7, ♦ 

<e POLLOWFO RY Y * 

rOM»iON/rC»^l/RX( 2 00) ,RY(2 00) , R 7 ( 200 ) , LCUT (2D0 ),LL,L0 
COMMON /P CR?/ RTF ,RnE2 ,R0E Y, AC ,PAC t PRTC, PTC, OFPTCt FPTCt 
2 FY.dfiv.TROT 

COMMON/C^ya/ix,! Y, I 7 , 1 X P , I Y P , I 7 -3 CX ,S CY , S C 7 , IDF FP 
C0'^MPN/crM7/Rl ,L?S,.SPX,SPZ,COX,COY,Cn7 
COMMON/CCyp/ I ox ,I OY, I D7 (200) , 0 XT , DYT , 07t ,T D ox , T POY, 

2 •^P(^7 

COMMON /COM2? /XLL ,XLF , 7LL ,71 F 

♦ TH'^ FOLIOWIK'D CCNSTANTc aRE USED oyTHEP RY THI<; SUR- ♦ 

♦ POUTINF TD GENFRATF thp LATTICE OR IN CTHFR SURROUT INFS * 

♦ AS RFOMTPFD . ♦ 



XLL=C. 6=7 
XLF = 2. 131 
7LI =C. =0 
7LF=3.Cl 
M = 2 

Tox= 1 
IDY=23 
I X=5 
I Y=5 
I7 = E 

= 0=1.C/50RT( 2.0) 

SC X=RO 
=C Y=RC) 

SC7 = 1 .0 
nXT = 2. C*=C X 
0YT = 2. C*SC Y 

07 T=? . c*FC 7 
OX=PP 
OY=0.0 
D7=0. E 

* THE FOLLOWING SERIFS QF NESTFn pg LOOPS GENFR&TF THE ♦ 

♦ lattice and insure IT CONTAINS THE REQUIRED SYMMgTPY * 

^ :C< 4c 4c <t 4 e ^ ’(t A 4 4c # <c * 4c « <1 <> X: A <> <c 4c ^ ^ 4< iCt 4 e * * 4c # 4> « ^ 4c 4‘ ^ 



)T = 0 
Y=-SC Y 

no 60 J=1,TY 
Y=Y+SC Y 
KT=C 



7=7 + SC7 
IT = 0 
Xc=- SC X 
MA =M 

DO AC T=l,IX 
X=X+ SC X 

IF( IT-( IT/2) *2) 

11 IF( JT-( JT/2) *2) 

12 TF(KT-(KT/2) *2) 

21 IF( .IT-( JT/2) *2) 

22 IF{KT-(kt/ 2) *2) 
30 P X( M) = x 

R Y( M)=Y 
R7( M)=7 
L=M 
M = M+ 1 

^0 IT=IT+1 
mb=m-i 



21 ,11 ,21 
40,12 ,40 
40 , 30,40 
22 , 40,22 
30 , 40,30 
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IF(mB-mA) 21C,21C,1<58 
I-*? JT.cq. 1 ) r,5 Tn 21'" 

IF(JT.F0.3) GO TO 21 B 
icq no 20C 1=^4 , mb 

ic( 3.9*Sr X-»X( ! n 200,200,202 
2C2 IF{ 7-P7( ! ) ) 200,200,203 
203 PX{ M)=PX{ I )4-0X 
RY( m)=fY( I )4-0Y 
R7(M)=R7( I )i-07 

M = M+1 

200 CO'MTINIIF 
GO TO 210 
21C lp=m-2 

RX( M ) =R x( LP ) -nx 

PY( M) =PY( LP) +DY 
P7 ( M) =P7 ( LP ) +D7 
M = M+ 1 
GO TP ICC 
?10 CnMTIN'UF 
CQ kt=kt+i 
jt=jt+i 

LL=m- 1 
00 2 I=2,LL 
2 I07(I)== 

♦ THF FPLL'’WTK:r, pr loop RTVIOFS ’’MF OPYSTAL INTO 4N ACTIVE* 

* PCPF NP AN' IN'CTTVP cyppouNPI mg "0 = 11. * 

P'^ c T = 2,LL 
LPUT( I >=C 

IF(PY( T ) .G-".2.2) LCnT{I)=l 

IF(P7( I ).GT.3.c) LCUT(I)=1 

TP(P7( 1 ) .LT.0.4) LCUTm=l 
c CO^lTI^'l!^ 

LCUT( n=f 

Rl=(Cnc+PY(2))/4Br(CrY) 

Rcryv V 
CMP 
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Computet Program - SUBROUTINE DLlll 



^cif*itc*^itft^********** SURROUTINF DLlll ♦♦♦♦♦♦♦♦♦♦♦**♦♦*♦♦♦♦♦ 

* THIS Tc; A IATTICE nENERATQR FDR THF DIAMOND LATT ICF IN ♦ 

* THE (111) orientation. ♦ 

* THF Crystal develops in the order x, followed by z, ♦ 

* FOLLOWED RY Y ♦ 

i/: 4:***-*f******^***** 4 ********^*************t******** ********** 

COMMON /CCM1/RX( 200) ,RY(2 00) . R Z ( 200 ) , LCUT ( 200 ) ,t.L,LD 
C0MMCN/CCM2/R0E ,R0E2 ,ROEM, AC ,PAr,PPTC,PTC,PFPTC, FPTC» 

2 pm^dptv.TPOT 

COMMON./cf ma/i x,I Y,I 7,1 XP,IYP, IZPfSCX ,SCY,SC7, IDEEP 
common /CCM 7/RI ,LSS ,SDX,SOZ ,COX,COY ,COZ 
COMMON/rrM^;/ IDX.IDY, 107(200) , DXT , DY T , 0 ZT , T PCX, T poy , 

2 TPP7 

common 7C0M22 7XLL ,XLF ,7LL ,7LF ,YLL,YLF 
************************************************************ 

* the FOLLOWING CONSTANTS ARE OPED EITHER BY THIS SUB- ♦ 

* ROUTINE TO generate THE LATTICE OR IN OTHER SUBROUTINES ♦ 

* AS REOUIRFO, ♦ 

************************************************************ 

XLL=0. 

XLF=2. 131 
YLL=2. 5 

YLF=i.4 
ZLL=1. 21 E 
ZLF=3.6F4 
M=2 
LA = 2 
IB = l 
IDX=l 
TDY=EE 

DO 2 1=2, lEC 
2 I07(T)=5 
IX=5 
IY=7 
TZ = E 

SCX=1.C/S0RT(2.0) 

SCY=2.C/S0RT(3.0) 

SC7=SQRT( l.S) 

SSC7 = SC7 /3.0 
nXT=2.C*SC X 
DYT=2.C*SORT(3.0) 

DZT=2.C*SC7 

DX=0.C 

OY=SORT( C. 75) 

07=0.0 

****************** 4 ***************************************** 

« the following series CIP NESTED DO LOOPS GENERATE THE * 

* lattice AN>^ insure it CONTAINS TH= REQUIRED SYMMETRY ♦ 

************************************************************ 

JT=0 
Y=- SC Y 

DO 50 J=l,IY 
Y=Y+SC Y 
KT = 0 

JTS = JT + .)T/3 
7=- CC7 

DO 50 K=1,I7 
IF(K.EO.E) GC to 80 
°0 CONTINUE 
Z = 7 + SC7 
IT = 0 
X=- SC X 

DO 40 1 = 1, TX 

X=X+SC X 

1N = IT<- JTS + KT 

IF( TN-f IN/2)>*2) 40,30,40 
30 R X( M) = x 
RY(M)=Y 
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TC( jt-3*( JT/3) ) 31,35,31 

31 JTT=JT 

32 JTT=JTT-3 
IF(JTT) 33,35,32 

33 JTT=JTT+3 
7P=JTT 

R7{M)=7+7P*$SCZ 

r,n m 37 

35 R7(M)=7 
37 L=M 
M=M+ 1 

no TO 4C 

PO IF( J.50. n 00 TO «0 

IF(.I.FO.<) GC TO QO 

TF( J.F0.7) r,C TO 90 

GO TO P5 
4C IT=IT+1 
cc kt=KT+1 

P5 IF( J.PQ.7) G': tq 100 
=0 

on 200 M=L«.,L 

IC( X-P X{ T.M *SCX) 2 00 ,2 00 ,202 

202 IF( a.C*5C7-p 7( IA n 200,203,203 
2C3 RX( L + IP )=RX( IA )+0X 
RY( L+IP )=RY( !A )+DY 
R7(L+IP)=PT( IA)+07 
IB=IP + 1 

MA=M.A + I 

20c COMTINUF 
y=L+WA+ 1 
LA = M 
IB=1 

fC JT=.)T+1 
ICC LL=M-i 

# THF FOLLOViTMO OC LOOP DIVIOFS th^ CRYSTAL INTO AN ACTIVE* 

* CORF ANO am TKA.CTIVF FURRGUNDTNG CHFLL. * 

4c :«c 4c i{c Itc !(ni( !(c X: ^ 4 : :4c ^ 4: ^ Xt 4: * 4 : Xc Xc # 4: ^ ^ * * # * « ^ * * Xc 4c Xc # XI ^ 4c Xc 

nn 3 I=?,lL 
LCUT( T )=C 

IF(PX(I).LT.C.7.rR.RXm.GT.2.2) LCUTI I ) = 1 
TF(P7m.LT.1.4.rR.R7(I).G’^.3.-^) LCUT( I )=1 
IF(P Y( I ) .GT. 5.2 ) LCUT(n=l 
3 CONTTNLF 
LCi'T( n=c 

P l=( pnc+p Y» 2 ) ) /ABS(CnY) 

P c T|jp (g 
FNO 
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Computer Program 



Exponential Interpolator 



^cm*l^************ EXPONENTIAL INT^RPCLATnp ♦♦♦♦♦♦♦♦*♦*♦♦♦♦♦ 

♦ READ LATTICF-LATTirE POTENTIALS AND DETERMINE 1= A BORN-* 

♦ MAYER ftPORT XI VA TION IS NECESSARY ♦ 

READ( 5 ,982S) ( ITI TLE ( I ) ,1=1 ,20) 

REAO( 5,c«lC) (TPNT(I) ,1=1 ,2 0) 

READ(S,cpiS)( TPNT( I ) , I =21 ,500) 

REAn(5,<^a2C)KA,TPCNST 
IP(KA) 6C5,4C5,410 
ACS ICHFrK=4AC 

GO 'p cif- 

« CnNSTPUr.T A BORN-mayER pCTENTIAL * 

410 DP 41" I=KA,500 

415 TPNT( I )=TPNT( KA ) *FXP(TPCNST*( ( I /lOO . )-(KA/ 100 . ) ) ) 

♦ 4c 4c 4i 4c :(( * 4t :([ <t ](( <( i(t 4i i(t * A 4< 4i 41 4< 4c 4c 4c 4< * # 4t 4 e <c 4ci(c 4c i(c ;(< i4c 4n(i 4i 4c « itc it( )(n4< 4c 4c 4c « 

♦ make an EXPCN’ENTIAL I nitFRPCLATI on between C^LCULATFC ♦ 

♦ Phttn^ialS * 

ICHFf K=KA 

416 ICnilNT = 2C 

417 DO 41C 1=1, 2C 
NCnUNT=irC(JNT+20 

BDFLX=-ALC0( TPNT(NCQI)NT) /TPNT(lcnUNT ) )/2C. 

419 TPNT( !CrUNT+I )=TPNT( IC0 UNT)c4cFXO ( (-1 ) *BDELX ) 

IF( ICHCCM-NCODNT) 422 ,422,421 

421 IC0UNT=NCPLN7 
r,p TP 4 17 

4c4c4c4c4c4<4c4c4c4c4c44c^c^444c:^4c4c4c4c4c4c4c4c4c4:4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c 

♦ DETERMINE POTENTIAL AT NEAREST NEIDHBOR OIS T A NC E , ? FRO * 

♦ POTENTIAL table, AND CONVERT POTENTIAL TO ELECTRON VOLTS ♦ 

4c4c4c4c4c4c4c4c4c4c44c4c4c44c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4<4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4(4c4c4c4c4c4c4c 

4 22 TTPNT = TPNT( I<=C ) + { ( TPNTI I FC+1 ) -T PNT ( I FC ) )* ( AFC-IFC ) ) 
no 420 1=1,500 
TPNT( I ):=( TDNT( I )-TTPNT) 4C13.6 
IF( TDNT( 1 ) )425,420,420 
425 TPNT(I)=C.O 

420 CONTINUE 

4c4'4c4c4c4c4c:^4c4cc^4<4c4c4c4c^4c4c4c4c4c4‘4>4c4c4cc^4c4cA4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4:4c4c4c4c 4c A* 4c 4c 4c 4c 4c 4c 

♦ CnNSTRlICT A LATTICE-LATTICE FORCE ^ABLE ♦ 

4c 4c 4c 4c 4c 4< 4c 4c 4c 4c 4c 4c 4c 4c 4< 4< 4c 4c 4> 4c 4c 4c 4 4< 4c 4c 4c « 4c 4c 4c 4c 4c 4c 4< 4c 4c 4< 4c 4c 4c 4c 4c 4c 4( 4< 4> 4c 4> 4c 4c 4c « 4c 4> 4c 4c 4c 4c 4c 

DO A30 1 = 1 ,<^99 

4 30 TEORC = ( I )=( TPNT(I ) -TPNT(l + l) ) 4cCVE0A 

4c4c4c4c4c4c4c4c4c4c4c4c4c4c^4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c 

♦ P^AP bullet-lattice POTENTIALS AND DETERMINE IF A BORN- ♦ 

4c maYFR APPROXI va TI on IS NECESSARY ♦ 

4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c « 4c 4c 4c « 4c 

R EA0(5,c 625 HI TITLE (I) ,1=21,40) 
RFAD(5,cpiC)(BPNTm,I=l,2 0) 

PEAD( 5,C81") (BPNTd ) ,1=21 ,500) 

RFAD( 5 ,ce2C) KA ,BPCNST 
I<=(KA) 4A5,445,450 
445 ICHECK=46C 
on TP 456 

4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c44c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 

♦ CONSTRUCT A BORN-mAYER POTENTIAL ♦ 

4c4c^^4c4c4c4c4c4c4c44c4c4c4c4c4c4c4c4c4c44c4c4c4c4>4c4c4c4c4c4c4c4c4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 

450 DP 455 I=KA,500 

45 5 BPMT( I )=BPNT(KA ) 4cE X° ( BPCNST4c ( ( I / 100 . ) - ( K A/ 100 . ) ) ) 

4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 

4c make an exponential INTERPOLATION BETWEEN CALCULATED * 

♦ potentials ♦ 

4c4c4c4c4t«4c4c4c4c4c4c4c4c4c4c4c4c4i4c4c4c4t4c4c4c4c4i4t4t4i4c4c4c4c4c4c4c4c4t4c4c4c4c4c4c4c4c4c4c4t4c4c4c*4c4c4c4c4c 

ICHFCK=KA 

456 ICPUNT=2C 

457 DO 459 1=1,20 
NC0UNT=TCCLNT+20 
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f^OFlX^-A I BPNT( NCnUNTJ /R PMT ( T CD IMT ) ) / 20 . 

RPNT( !Cri'WT+ I ) =ROmt( I COUNT ) *EXC ( (- I) ♦BDFLX ) 

TF( TCHPCK-NCCUNT) 462 ,462,461 

461 TrnuNT=Nr rijNT 

r.n TP AC 7 

* DFtfdmtm= PPTc,jt^^l at NFARF'^T NFIOHROO DI ST ANC F , 7 ERO ♦ 

* POT=MTTyfL ■^4RL=,AND convert PrTENTIAL TC ELFC^PQN VOLTS ♦ 

462 RRPNT^BCMTf I FC » + ( ( B ^N T ( I FC + 1 ) -B P\'T ( I FC ) ) * ( 4F C- T FC) ) 
nn 660 1 = 1, SCO 

rpnt(I)=(RPnT(I)-RRPMT)«13.6 
IF(RPMT( I ) >475,460,480 
67C RPNT( I )=C.C 
4?C CCN'^INUF 

« <i 4c 4 « :4c 4i 4c ^ 4c ^ 4 <c 4 « 4c ^ 4 A 4c 4c :(( 4c ;t( 4c A 4c ^ # A <C ^ 

4c CONSTRUCT A BULLET-L ATTI CF FORCE TABLE * 

4c:^4c4c4c4c4c4c4c4c4c4<4c4c4444i44c4c4c4c4ci^4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4(4c>^4c4<4c4c4c4c4c4c4c4c4c4c4c4c4c 
4,65 on 6C0 Ic=1,4<?Q 

4CC BFpPCE ( I ) = (RCNT( I )-RP\T( 1 + 1 ) ) 4<r VEHA 
4c4c*4c*4c4c4c4c4c4c4c4c4c4t4c4c END OF I NT E R P F>1 A T 0 ® 4: * 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4c 4< 4c 4c 4c 4c 4c 4c 
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SILICON (100) SURFACE 






Figur e I . 
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S ILICON (110) SUFACE 






Figure 2, 
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SILIC0N<I1I> SURFACE 




Figure 3 * 
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Figure 4. The shape and sizes of the major 
channels in silicon. 
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RESIDUAL ACTIVITY (%) 




Figure 5* Range distribution curves for 40 KeV 
125 ,, . . , - 

Xe in silicon. 

Note: p4.9/cm^ = k2% = 15.5LU. 

•{Ref 1, Fig. 5) 
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<lOO>CHANNEL 
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POTENTIAL ENERGY (RYDBERGS) 



Xe-Sl POTENTIAL FUNCTIONS 




Figure 10. 
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Figure 12. 
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Figure l4 • 
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Figure l6 . 
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Figure 17 - 
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Figure 21. 
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SILICON<lh CHANNEL 




Figure 22. 
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Figure 23. 
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Figure 26. 
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